library(sf)
## Warning: package 'sf' was built under R version 4.3.2
## Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE
library(rsm)
## Warning: package 'rsm' was built under R version 4.3.1
library(broom)
## Warning: package 'broom' was built under R version 4.3.1
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.2
options(scipen = 1000)

Read data

base_dir = getwd()

data_dir = paste0(base_dir, "/", "data_maplayers_one_field.csv")

data_wkt <- st_read(data_dir, options = "GEOM_POSSIBLE_NAMES=WKT")
## options:        GEOM_POSSIBLE_NAMES=WKT 
## Reading layer `data_maplayers_one_field' from data source 
##   `C:\Users\BX2A428\JD\Projects\R\In-season-EONR-estimation-with-smart-farming-data\data_maplayers_one_field.csv' 
##   using driver `CSV'
## Simple feature collection with 1840 features and 66 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 430277.1 ymin: 5501739 xmax: 430537.1 ymax: 5502069
## CRS:           NA
maplayers_sf <- data_wkt[,!(names(data_wkt) == "WKT")] 
st_crs(maplayers_sf) <- st_crs("+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs")


## Covert datatype to numeric
maplayers_sf[, sapply(maplayers_sf, is.character)] <- lapply(maplayers_sf[, sapply(maplayers_sf, is.character)], as.numeric)
## Warning in `[<-.data.frame`(`*tmp*`, , sapply(maplayers_sf, is.character), :
## provided 66 variables to replace 65 variables
maplayers_df = sf::st_drop_geometry(maplayers_sf)

Yield-applied N relationship using Quadratic function

#### Wet-Mass vs applied_rate, N-treatment
ggplot(maplayers_df, aes(x = applied_rate, y=wet_mass_idw)) + geom_point(size=3)  +
                      geom_smooth(method = "lm", formula = y ~ poly(x,2), se=FALSE) + 
                      ggtitle("quadratic regression") + xlab("N-applied, kg/ha") + ylab("harvested grain, kg/ha")

Filter the best identified multispectral indicess during feature selection process using ANOVA and Lasso regression.

candidate_VIs = c(
                'cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632',
                'cloud.coverage.0_LCI_2020.03.25.10.36.59_7632',
                'cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632',
                'cloud.coverage.0_CRI550_2020.03.25.10.36.59_7632',
                'cloud.coverage.0_EVI_2020.03.25.10.36.59_7632',
                
                 'cloud.coverage.0_CARI_2020.03.30.10.37.00_8331',
                'cloud.coverage.0_CARI2_2020.03.30.10.37.00_8331',
                'cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331',
                'cloud.coverage.0_NDRE_2020.03.30.10.37.00_8331',
                
                'cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536',
                'cloud.coverage.0_LCI_2020.04.04.10.36.59_1536',
                'cloud.coverage.0_CRI550_2020.04.04.10.36.59_1536',
                'cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536',
                'cloud.coverage.0_NDRE_2020.04.04.10.36.59_1536',
                "applied_rate",
                "wet_mass_idw")




candidate_VIs_df = subset(maplayers_df, select = candidate_VIs)

Modeling with RSM

reduced_model_formula = wet_mass_idw ~ SO(applied_rate,
                        cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632,
                        cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632,
                        
                        cloud.coverage.0_CARI_2020.03.30.10.37.00_8331,
                        cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331,
                        
                        cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536,
                        cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536)




full_model <- rsm(reduced_model_formula, data = candidate_VIs_df)
mdl_summary = summary(full_model)
mdl_summary
## 
## Call:
## rsm(formula = wet_mass_idw ~ SO(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, 
##     cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, 
##     cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, 
##     cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536), data = candidate_VIs_df)
## 
##                                                                                                         Estimate
## (Intercept)                                                                                          -8665.28562
## applied_rate                                                                                           277.92723
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632                                                     327883.37382
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632                                                      -52992.87645
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331                                                       48258.42439
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331                                                      1813.89118
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536                                                    -235311.44334
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536                                                      5941.49507
## applied_rate:cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632                                          6156.23351
## applied_rate:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632                                             72.03097
## applied_rate:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331                                            704.62688
## applied_rate:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331                                           11.23289
## applied_rate:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536                                         -7282.89639
## applied_rate:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536                                           21.91666
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632      362872.49773
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331     -583776.08334
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331     93856.88086
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536   -7013590.66180
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536   -303837.76290
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331        -7744.49312
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331       -910.07235
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536       51900.43518
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536       3275.02536
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331      15020.05939
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536    -3492702.94031
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536      12528.83245
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536   -143638.55208
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536      854.88002
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536    200675.67714
## applied_rate^2                                                                                          -2.55804
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632^2                                                  2697229.06188
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632^2                                                    -16024.82153
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331^2                                                   1023909.60102
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331^2                                                    -671.31560
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536^2                                                  7523218.75660
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536^2                                                     721.81252
##                                                                                                       Std. Error
## (Intercept)                                                                                           9348.28152
## applied_rate                                                                                           156.82393
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632                                                     292493.87198
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632                                                       13314.31653
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331                                                       60808.30099
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331                                                      2575.15960
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536                                                     310413.34477
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536                                                      3552.04074
## applied_rate:cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632                                          2250.78720
## applied_rate:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632                                             93.33038
## applied_rate:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331                                            606.04170
## applied_rate:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331                                           19.62163
## applied_rate:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536                                          2322.32627
## applied_rate:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536                                           27.74469
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632      172054.07224
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331     1102246.45955
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331     42395.01989
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536    6038832.45311
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536     65441.75720
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331        43744.44000
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331       1500.73728
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536      194774.54559
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536       2532.93673
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331      11098.55669
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536     1352880.15990
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536      12301.17326
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536     47821.72872
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536      542.65633
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536     63366.54314
## applied_rate^2                                                                                           0.86505
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632^2                                                  2691649.96891
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632^2                                                      3573.73923
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331^2                                                    153496.48236
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331^2                                                     287.43781
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536^2                                                  4018842.95373
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536^2                                                     496.64474
##                                                                                                   t value
## (Intercept)                                                                                       -0.9269
## applied_rate                                                                                       1.7722
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632                                                    1.1210
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632                                                    -3.9801
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331                                                     0.7936
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331                                                   0.7044
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536                                                   -0.7581
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536                                                   1.6727
## applied_rate:cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632                                       2.7351
## applied_rate:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632                                        0.7718
## applied_rate:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331                                        1.1627
## applied_rate:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331                                      0.5725
## applied_rate:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536                                      -3.1360
## applied_rate:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536                                      0.7899
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632     2.1091
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331    -0.5296
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331   2.2139
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536   -1.1614
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536  -4.6429
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331     -0.1770
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331   -0.6064
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536     0.2665
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536    1.2930
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331    1.3533
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536    -2.5817
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536    1.0185
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536  -3.0036
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536  1.5754
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536   3.1669
## applied_rate^2                                                                                    -2.9571
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632^2                                                  1.0021
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632^2                                                  -4.4840
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331^2                                                   6.6706
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331^2                                                -2.3355
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536^2                                                  1.8720
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536^2                                                 1.4534
##                                                                                                           Pr(>|t|)
## (Intercept)                                                                                               0.354082
## applied_rate                                                                                              0.076526
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632                                                           0.262440
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632                                                    0.00007160280313
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331                                                            0.427524
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331                                                          0.481287
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536                                                           0.448515
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536                                                          0.094560
## applied_rate:cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632                                              0.006296
## applied_rate:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632                                               0.440343
## applied_rate:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331                                               0.245117
## applied_rate:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331                                             0.567072
## applied_rate:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536                                              0.001740
## applied_rate:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536                                             0.429666
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632            0.035077
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331            0.596438
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331          0.026963
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536           0.245627
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536  0.00000368410010
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331             0.859497
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331           0.544314
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536            0.789912
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536           0.196185
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331           0.176118
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536            0.009910
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536           0.308574
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536          0.002705
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536         0.115348
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536          0.001566
## applied_rate^2                                                                                            0.003146
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632^2                                                         0.316443
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632^2                                                  0.00000778509314
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331^2                                                  0.00000000003378
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331^2                                                        0.019626
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536^2                                                         0.061370
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536^2                                                        0.146293
##                                                                                                      
## (Intercept)                                                                                          
## applied_rate                                                                                      .  
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632                                                      
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632                                                    ***
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331                                                       
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331                                                     
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536                                                      
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536                                                  .  
## applied_rate:cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632                                      ** 
## applied_rate:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632                                          
## applied_rate:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331                                          
## applied_rate:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331                                        
## applied_rate:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536                                      ** 
## applied_rate:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536                                        
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632    *  
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331       
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331  *  
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536      
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536  ***
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331        
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331      
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536       
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536      
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331      
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536    ** 
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536      
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536  ** 
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536    
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536  ** 
## applied_rate^2                                                                                    ** 
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632^2                                                    
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632^2                                                  ***
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331^2                                                  ***
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331^2                                                *  
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536^2                                                 .  
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536^2                                                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Multiple R-squared:  0.4066, Adjusted R-squared:  0.3951 
## F-statistic: 35.32 on 35 and 1804 DF,  p-value: < 0.00000000000000022
## 
## Analysis of Variance Table
## 
## Response: wet_mass_idw
##                                                                                                                                                                                                                                                                                                                           Df
## FO(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536)     7
## TWI(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536)   21
## PQ(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536)     7
## Residuals                                                                                                                                                                                                                                                                                                               1804
## Lack of fit                                                                                                                                                                                                                                                                                                             1804
## Pure error                                                                                                                                                                                                                                                                                                                 0
##                                                                                                                                                                                                                                                                                                                             Sum Sq
## FO(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536)  1024415381
## TWI(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536)  530111687
## PQ(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536)   158775478
## Residuals                                                                                                                                                                                                                                                                                                               2499956661
## Lack of fit                                                                                                                                                                                                                                                                                                             2499956661
## Pure error                                                                                                                                                                                                                                                                                                                       0
##                                                                                                                                                                                                                                                                                                                           Mean Sq
## FO(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536)  146345054
## TWI(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536)  25243414
## PQ(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536)   22682211
## Residuals                                                                                                                                                                                                                                                                                                                 1385785
## Lack of fit                                                                                                                                                                                                                                                                                                               1385785
## Pure error                                                                                                                                                                                                                                                                                                                    NaN
##                                                                                                                                                                                                                                                                                                                         F value
## FO(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536)  105.604
## TWI(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536)  18.216
## PQ(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536)   16.368
## Residuals                                                                                                                                                                                                                                                                                                                      
## Lack of fit                                                                                                                                                                                                                                                                                                                 NaN
## Pure error                                                                                                                                                                                                                                                                                                                     
##                                                                                                                                                                                                                                                                                                                                        Pr(>F)
## FO(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536)  < 0.00000000000000022
## TWI(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536) < 0.00000000000000022
## PQ(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536)  < 0.00000000000000022
## Residuals                                                                                                                                                                                                                                                                                                                                    
## Lack of fit                                                                                                                                                                                                                                                                                                                               NaN
## Pure error                                                                                                                                                                                                                                                                                                                                   
## 
## Stationary point of response surface:
##                                     applied_rate 
##                                       68.9526134 
##  cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632 
##                                        0.1404002 
##   cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632 
##                                        0.4542571 
##   cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 
##                                        0.2210777 
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 
##                                        0.1164771 
##  cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 
##                                        0.1401466 
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 
##                                        1.8876305 
## 
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 9605894.552479 1891690.024961   24357.837308       1.302358    -363.129399
## [6]   -3736.778354 -289463.272496
## 
## $vectors
##                                                           [,1]          [,2]
## applied_rate                                     -0.0004817841  0.0005765301
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632  -0.4412944672  0.6937323329
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632   -0.0058733357  0.0703332995
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331   -0.1644532894 -0.6800575429
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 -0.0088768487  0.0060642931
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536   0.8819535609  0.2216808467
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536  0.0160841100 -0.0461617788
##                                                          [,3]         [,4]
## applied_rate                                      0.004082907 0.9947852226
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632   0.177160547 0.0013894104
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632    0.856472383 0.0250415804
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331    0.342167820 0.0004836887
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 -0.053641044 0.0630152248
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536   0.163038067 0.0007407485
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 -0.297514600 0.0761685534
##                                                          [,5]        [,6]
## applied_rate                                     -0.098126742 -0.02745010
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632   0.015727534  0.01007352
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632    0.228642726  0.24255836
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331    0.013208515 -0.02810980
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331  0.802205000 -0.59110454
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536   0.009986534 -0.01855557
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536  0.542253845  0.76796576
##                                                          [,7]
## applied_rate                                      0.001661553
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632  -0.540608937
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632    0.386955620
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331   -0.626444871
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331  0.009935057
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536  -0.382085748
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 -0.139582105

Regression diagnostic plots

augmented_data <- augment(full_model)
## Warning: The `augment()` method for objects of class `rsm` is not maintained by the broom team, and is only supported through the `lm` tidier method. Please be cautious in interpreting and reporting broom output.
## 
## This warning is displayed once per session.
# Residuals vs Fitted plot
ggplot(data = augmented_data, aes(.fitted, .resid)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  ggtitle("Residuals vs Fitted") +
  xlab("Fitted Values") +
  ylab("Residuals") + labs(title = "Residuals vs Fitted Plot - Best-model")

# Normal Q-Q plot of residuals
ggplot(data = augmented_data, aes(sample = .std.resid)) +
  geom_qq() +
  geom_qq_line() +
  ggtitle("Normal Q-Q Plot of Residuals") + labs(title = "Normal Q-Q plot of residuals - Best-model")

## Compute model accuracy. RMSE percentage

# Calculate the squared differences between actual and predicted values
squared_errors <- (candidate_VIs_df$wet_mass_idw - full_model$fitted.values)^2

# Calculate the mean of squared errors
mean_squared_error <- mean(squared_errors)

# Calculate the RMSE (take the square root of the mean squared error)
rmse <- sqrt(mean_squared_error)

(rmse/mean(candidate_VIs_df$wet_mass_idw)) * 100
## [1] 14.58546

Compute N-opt

compute_n_opt = function(trainmodel, first.date.CARI2, first.date.ARI2, second.date.CARI, second.date.CRI550, third.date.CARI2, third.date.CRI700){
  
  
  mdl_summary = summary(trainmodel)
  # collected_terms = mdl_summary$coefficients[c(1:14, 30:36), ]
  
#######################################################################
    ### Linear terms plus estimated coefficient for N
    # c.0 = mdl_summary$coefficients[c(1, 3:8),]
#######################################################################  
  
  first_date_CARI2 = mdl_summary$coefficients[3] * first.date.CARI2
  
  first_date_ARI2 = mdl_summary$coefficients[4] * first.date.ARI2
  
  second_date_CARI = mdl_summary$coefficients[5] * second.date.CARI
  
  second_date_CRI550 = mdl_summary$coefficients[6] * second.date.CRI550
  
  third_date_CARI2 = mdl_summary$coefficients[7] * third.date.CARI2

  third_date_CRI700 = mdl_summary$coefficients[8] * third.date.CRI700

  c.0 =  mdl_summary$coefficients[1] + first_date_CARI2 + first_date_ARI2 + second_date_CARI + second_date_CRI550 + third_date_CARI2 + third_date_CRI700
  
  

#######################################################################
    ### Interaction terms plus estimated coefficient for N
    # C.1 = mdl_summary$coefficients[c(2, 9:14)]
#######################################################################  
  
  first_date_CARI2_interaction = mdl_summary$coefficients[9] * first.date.CARI2
  
  first_date_ARI2_interaction = mdl_summary$coefficients[10] * first.date.ARI2
  
  second_date_CARI_interaction = mdl_summary$coefficients[11] * second.date.CARI
  
  second_date_CRI550_interaction = mdl_summary$coefficients[12] * second.date.CRI550

  
  third_date_CARI2_interaction = mdl_summary$coefficients[13] * third.date.CARI2

  third_date_CRI700_interaction = mdl_summary$coefficients[14] * third.date.CRI700
  
  c.1 =  mdl_summary$coefficients[2] +  first_date_CARI2_interaction + first_date_ARI2_interaction + second_date_CARI_interaction + second_date_CRI550_interaction + third_date_CARI2_interaction + third_date_CRI700_interaction
  

  c.2 =  mdl_summary$coefficients[30]


  n_opt = (0.417/0.318) * (1/(2*c.2)) - (c.1/ (2*c.2))
  
  
  
  if(class(n_opt)== "data.frame"){
    
    
    names(n_opt) = "optimum_nitrogen"
    names(c.0) = "intercept"
    names(c.1) = "interaction_terms"

    n_opt$intercept = c.0$intercept
    n_opt$interaction = c.1$interaction_terms

    
    }
  
  
  return(n_opt)
  
  
  
}










#####################################################################################
################### Prediction of EONR ############################################
#####################################################################################



final_vis_predictors = c("cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632", 
                        "cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632", 
                         "cloud.coverage.0_CARI_2020.03.30.10.37.00_8331",
                        "cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331", 
                        "cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536", 
                        "cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536")


VIs_predictors_df = candidate_VIs_df[, final_vis_predictors]


n_inpute = compute_n_opt(trainmodel = full_model, 
             first.date.CARI2 = VIs_predictors_df[1], first.date.ARI2 = VIs_predictors_df[2],
             second.date.CARI = VIs_predictors_df[3], second.date.CRI550 = VIs_predictors_df[4],
             third.date.CARI2 = VIs_predictors_df[5], third.date.CRI700 = VIs_predictors_df[6])


hist(n_inpute$optimum_nitrogen,  main = "Economic optimum N rate", xlab = "N-kg/ha")

hist(candidate_VIs_df$applied_rate)

sf_final_data <- subset(maplayers_sf, select = final_vis_predictors)

sf_final_data$optimum_n = n_inpute[, 'optimum_nitrogen']

# Define a color palette with red, yellow, and green
color_palette <- colorRampPalette(c("red", "yellow", "green"))

# Create map view with continuous color scale
mapview::mapview(sf_final_data, zcol = "optimum_n", col.regions = color_palette,
        legend = TRUE, legend.style = "classic")

Calculate Total N

sf_final_data$optimum_n_total_N = sf_final_data$optimum_n  / 0.26


# Define a color palette with red, yellow, and green
color_palette <- colorRampPalette(c("red", "yellow", "green"))

# Create map view with continuous color scale
mapview::mapview(sf_final_data, zcol = "optimum_n_total_N", col.regions = color_palette,
        legend = TRUE, legend.style = "classic")

Plot two points with different EONR

The location indexes we used for visualization purpose are as follows, 1 and 1800

c_value_0 <- n_inpute$intercept[1]
c_value_1 <-n_inpute$interaction[1]
c_value_2 <- mdl_summary$coefficients[30]

# Single value N
N <- seq(30, 90, length.out = 100)

# Quadratic function
Y <- c_value_0 + c_value_1 * N + c_value_2 * N^2


# Calculate N for the maximum
N_max <- -c_value_1 / (2 * c_value_2)

# Calculate corresponding f(N) for the maximum
Y_max <- c_value_0 + c_value_1 * N_max + c_value_2 * N_max^2

# Plot the curve
plot(N, Y, type = "l", col = "blue", xlab = "N kg/ha", ylab = "Y(N) kg/ha", main = "Quadratic Function")

# Add a vertical line at the maximum
abline(v = N_max, col = "red", lty = 2)

# Add a legend
legend("bottomright", legend = expression(Y(N) == c[0] + c[1] * N + c[2] * N^2), col = "blue", lty = 1)

# Add text annotations for c_value_0, c_value_1, and c_value_2
text(N[1], Y[1], paste("c[0] =", round(c_value_0, 2)), pos = 4, col = "black")
text(N[50], Y[50], paste("c[1] =", round(c_value_1, 2)), pos = 2, col = "black")
text(N[80], Y[80], paste("c[2] =", round(c_value_2, 2)), pos = 1, col = "black")

# Print the results
cat("N for maximum:", N_max, "\n")
## N for maximum: 54.48445
cat("Maximum value Y(N):", Y_max, "\n")
## Maximum value Y(N): 5804.641
c_value_0
## [1] -1789.041
c_value_1
## [1] 278.7468
c_value_2
## [1] -2.55804
c_value_0 <- n_inpute$intercept[1800]
c_value_1 <-n_inpute$interaction[1800]
c_value_2 <- mdl_summary$coefficients[30]

# Single value N
N <- seq(30, 90, length.out = 100)

# Quadratic function
Y <- c_value_0 + c_value_1 * N + c_value_2 * N^2


# Calculate N for the maximum
N_max <- -c_value_1 / (2 * c_value_2)

# Calculate corresponding f(N) for the maximum
Y_max <- c_value_0 + c_value_1 * N_max + c_value_2 * N_max^2

# Plot the curve
plot(N, Y, type = "l", col = "blue", xlab = "N kg/ha", ylab = "Y(N) kg/ha", main = "Quadratic Function")

# Add a vertical line at the maximum
abline(v = N_max, col = "red", lty = 2)

# Add a legend
legend("bottomright", legend = expression(Y(N) == c[0] + c[1] * N + c[2] * N^2), col = "blue", lty = 1)

# Add text annotations for c_value_0, c_value_1, and c_value_2
text(N[1], Y[1], paste("c[0] =", round(c_value_0, 2)), pos = 4, col = "black")
text(N[50], Y[50], paste("c[1] =", round(c_value_1, 2)), pos = 2, col = "black")
text(N[80], Y[80], paste("c[2] =", round(c_value_2, 2)), pos = 1, col = "black")

# Print the results
cat("N for maximum:", N_max, "\n")
## N for maximum: 69.87685
cat("Maximum value Y(N):", Y_max, "\n")
## Maximum value Y(N): 10389.38
c_value_0
## [1] -2100.952
c_value_1
## [1] 357.4955
c_value_2
## [1] -2.55804

Plot a bundle of points.

Picking a few locations in the field and ploting the N-yield response.

Candidate locations are 1, 1800, 100, 200, 56, and 440.

c_0 <- n_inpute$intercept
c_1 <-n_inpute$interaction
c_2 <- mdl_summary$coefficients[30]

c_0_subset = c_0[c(1, 1800, 100, 200, 56, 440)]
c_1_subset = c_1[c(1, 1800, 100, 200,  56, 440)]

c_values_0 <- c_0_subset
c_values_1 <-c_1_subset
c_value_2 <- c_2

# Single value N
N <- seq(0, 100, length.out = 400)

# Create an empty plot with specified ylim
plot(N, rep(NA, length(N)), type = "n", col = "white", xlab = "N kg/ha", ylab = "Y(N) kg/ha", main = "Quadratic Functions", ylim = c(1000, 12000))

# Plot the curves
for (i in 1:length(c_values_0)) {
  c_value_0 <- c_values_0[i]
  c_value_1 <- c_values_1[i]
  
  # Quadratic function
  f <- c_value_0 + c_value_1 * N + c_value_2 * N^2
  
  # Add the curve to the plot
  lines(N, f, col = i + 1)  # Use different colors for each curve
}

# Add a legend
legend("topright", legend = c(expression(Y(N) == c[0] + c[1] * N + c[2] * N^2), paste("Set ", 2:5)), col = 2:(length(c_values_0) + 1), lty = 1)